library("knitr")
opts_chunk$set(
# cache=FALSE,
message=FALSE,warning=FALSE)
library("ggplot2") # для построения графиков
library("rasterVis")
library("fields")
library("deSolve")
library("bvpSolve")
Пакет rasterVis предназначен для изображения данных на реальных географических картах, поэтому там нужно понятие проекции. Мы пока просто введем это шаманское заклинание
proj <- CRS('+proj=longlat +datum=WGS84')
\[ \left\{ \begin{array}{l} \dot{y}_1=y_2 \\ \dot{y}_2=y_1+\cos(y_2) \end{array} \right. \]
Задаем решетку, и рассчитываем \(\dot{y}_1\) и \(\dot{y}_2\) в точках решетки:
y1 <- seq(-6, 6, .05)
y2 <- seq(-6, 6, .05)
df <- expand.grid(y1=y1, y2=y2)
df$y1dot <- df$y2
df$y2dot <- df$y1+cos(df$y2)
Рассчитываем длины и углы для стрелочек, помещаем результат в объект Raster.
df$len <- sqrt(df$y1dot^2+df$y2dot^2)
df$angle <- atan2(df$y1dot,df$y2dot)
df2 <- df[c("y1","y2","len","angle")]
rast <- rasterFromXYZ(df2,crs=proj)
Строим классический график со стрелочками
vectorplot(rast,isField=TRUE)
Строим няку с капельками
streamplot(rast, isField=TRUE)
Простой график можно руками построить без доп. пакетов. При этом нам нужно самостоятельно уменьшить количество стрелочек.
y1 <- seq(-6, 6, .5)
y2 <- seq(-6, 6, .5)
df <- expand.grid(y1=y1, y2=y2)
df$y1dot <- df$y2
df$y2dot <- df$y1+cos(df$y2)
plot(df$y1,df$y2,pch=".")
arrow.plot( df$y1,df$y2,df$y1dot,df$y2dot,
arrow.ex=0.03,length=0.05)
Описываем саму систему:
eq1 <- function(t,y,parampampam) {
return(list(c(
y[2],
y[1]+cos(y[2])
)))
}
Начальные условия:
y.start <- c(y1=1,y2=4)
Точки, в которых компьютер будет считать функцию:
t <- seq(0,10,by=0.01)
Решаем
sol <- ode(y=y.start,times=t,func=eq1)
sol <- data.frame(sol)
head(sol)
## time y1 y2
## 1 0.00 1.000000 4.000000
## 2 0.01 1.040018 4.003678
## 3 0.02 1.080076 4.007785
## 4 0.03 1.120176 4.012326
## 5 0.04 1.160324 4.017305
## 6 0.05 1.200524 4.022725
qplot(data=sol,time,y1)
Функция ode возвращает матрицу, а для рисования графиков удобнее табличка с данными, data.frame. Строчка sol <- data.frame(sol) переделывает матрицу в таблицу с данными.
Описываем саму систему:
eq1 <- function(t,y,parampampam) {
return(list(c(
y[2],
y[1]+cos(y[2])
)))
}
Граничные условия:
y.start <- c(y1=1,y2=NA)
y.final <- c(y1=42,y2=NA)
Точки, в которых компьютер будет считать функцию:
t <- seq(0,10,by=0.01)
Решаем
sol <- bvptwp(yini=y.start,yend=y.final,
x=t,func=eq1,
nmax=2000)
sol <- data.frame(sol)
head(sol)
## x y1 y2
## 1 0.00 1.0000000 -1.553150
## 2 0.01 0.9845193 -1.543001
## 3 0.02 0.9691398 -1.532904
## 4 0.03 0.9538610 -1.522860
## 5 0.04 0.9386824 -1.512868
## 6 0.05 0.9236035 -1.502928
qplot(data=sol,x,y1)
Есть несколько способов представить себе функцию от двух переменных, \(z(x,y)\):
Создаем data.frame с декартовым произведением двух векторов
df <- expand.grid(x=seq(-2, 2, .01), y=seq(-2, 2, .01))
Изобразим функцию \(z(x,y)=(3\cdot x^2+y)\cdot e^{-x^2-y^2}\).
Cоздаем переменную z как функцию от x и y
df$z <- with((3*x^2 + y)*exp(-x^2-y^2),data=df)
r <- rasterFromXYZ(df, crs=proj)
Линии уровня функции z
contour(r)
Капельки текущие по градиенту
streamplot(r)
Направление градиентов, заодно вид сбоку для графика функции
vectorplot(r)